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Abstract 


The  advantages  of  solving  potential  problems  using  an  overdetermined  boundary 
integral  element  method  are  examined.  Representing  a  two-dimensional  potential 
solution  by  an  analytic  complex  function  forms  two  algebraic  systems  from  the  real 
and  imaginary  parts  of  the  discretized  form  of  Cauchy’s  theorem.  Depending  on 
which  boundary  condition  is  prescribed,  the  real  or  the  imaginary  algebraic  system 
is  diagonally  dominant.  Computations  show  that  the  errors  of  the  strong  system 
(diagonally  dominant)  often  have  almost  the  same  value  as  those  of  weak  system 
(diagonally  non-dominant)  but  with  the  opposite  sign.  The  overdetermined  sys¬ 
tem.  composed  of  the  combination  of  the  real  and  imaginary  parts,  tends  to  average 
these  errors,  especially  for  circular  contours.  An  error  analysis  and  convergence 
studies  for  several  geometries  and  boundary  conditions  are  performed.  A  method¬ 
ology  for  handling  computational  difficulties  with  contour  comers  is  outlined.  Fur¬ 
ther  modifications  are  proposed  and  tested  that  show  exponential  convergence  for 
smooth  contours. 
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1  Introduction 


The  boundary  integral  method  (also  called  the  panel  method)  is  powerful  for  solving  potential 
problems.  In  this  method,  the  computational  domain  becomes  the  enclosing  boundary,  reducing 
the  effective  dimension  of  the  problem  by  one.  This  significantly  decreases  the  computational 
effort  unless  the  contour  is  highly  contorted. 

Boundary  integral  methods  are  usually  derived  from  Green’s  theorem  with  an  appropriate 
free-space  Green’s  function.  However,  the  solutions  of  two-dimensional  potential  problems  dis¬ 
cussed  here  can  be  described  by  a  complex  analytic  function.  This  allows  the  boundary  integral 
method  to  be  formulated  from  Cauchy’s  integral  theorem  (Vinje  and  Brevig  [1],  Lai  and  Hro- 
madka  [2]).  The  elegance  and  simplicity  of  complex  analysis  carries  over  to  computations  as  well. 
While  the  Green’s  function  formulation  and  the  complex  method  are  not  directly  compared  here, 
Dold  and  Peregrine  [3],  indicate  that  the  latter  method  (from  [1])  is  clearly  superior.  Schultz  [4] 
shows  that  the  error  of  the  complex  formulation  can  be  further  reduced  by  solving  the  resulting 
two  systems  of  algebraic  equations  in  a  least-squares  sense  for  nonlinear  breaking  wave  simula¬ 
tions.  In  this  paper,  we  expand  and  formalize  these  ideas  for  model  potential  problems. 

In  general,  an  integral  equation  is  solved  by  a  numerical  model  that  assumes  the  boundary  is 
composed  of  piecewise-polynomial  curves  (panels)  and  the  known  and  unknown  boundary  values 
ue  approximated  as  piecewise-continuous  functions  along  the  boundary.  The  truncation  error  in 
solving  an  integral  equation  is  determined  by  the  approximations  of  the  boundary  shape  and  the 
known  and  unknown  function  values  on  the  boundary.  A  rather  extensive  error  analysis  of  the 
conventional  boundary  integral  method  for  a  two-dimensional  Neumann  problem  is  given  by  Hess 
[5].  He  shows  that  using  parabolically-shaped  elements  with  linearly  varying  singularity  improves 
the  accuracy  when  the  effects  of  local  element  curvature  are  greater  than  that  of  the  function 
derivatives.  Unlike  the  Green’s  function  formulations,  the  complex  method  does  not  depend 
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on  the  shape  of  the  contour  between  nodes.  Most  complex  variable  boundary  integral  methods 
adopt  piecewise-linear  representations  of  the  complex  functions.  This  linear  interpolation  results 
in  second-order  accurate  integration  and  typically  gives  second-order  accurate  solutions  for  the 
boundary  integral  solution  as  well. 

The  accuracy  of  the  numerical  solution  is  affected  by  the  contour  shapes,  boundary  conditions 
and  nodal  spacings  as  well  as  the  possible  singular  behavior  of  the  resulting  solutions.  Moran  [6] 
shows  poor  convergence  for  a  sharp-angled  contour  (such  as  airfoils)  even  when  the  solution  there 
is  smooth.  In  addition,  singularities  usually  occur  at  an  abrupt  change  in  boundary  conditions  or 
at  boundary  corners  [12].  The  solution  accuracy  then  further  degrades  unless  special  precautions 
are  taken  such  as  supplementing  the  piecewise-linear  representations  with  singular  functions  [7]. 

In  this  paper,  the  error  is  analyzed  for  the  complex  variable  boundary  integral  method  using 
elements  with  linearly  varying  functions.  The  advantages  of  solving  both  sets  (real  and  imagi¬ 
nary)  of  algebraic  equations  together  in  a  least-squares  sense  are  discussed.  The  integration  error 
analysis  shows  the  correlation  between  solution  error  and  factors  such  as  geometric  curvature, 
nodal  spacing,  and  local  solution  gradients.  We  also  present  further  modifications  in  sections  6 
and  7  that  show  exponential  convergence  for  smooth  contours. 
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2  Cauchy  Boundary  Integral  Formulation 


Cauchy’s  theorem  is  used  to  derive  an  integral  equation  in  terms  of  the  complex  potential. 
The  complex  potential,  0,  is  analytic  inside  the  computational  domain,  R,  and  is  given  by 

0(z)  =  <i>+  itl>  ,  (1) 

where  the  two-dimensional  spatial  representation  is  given  by  z  =  x  iy,  and  ^  and  are  real 
functions  that  can  be  identified  as  the  velocity  potential  and  stream  function,  respectively,  for 
ideal  fluid  flow  problems. 

Figure  I  shows  the  problem  domain  and  bounding  surfaces.  Cauchy’s  theorem  gives 

/  =  «a/!?(Ct)  .  (2) 

JdR  «  ~  St 

where  a  is  0  or  2x  if  the  location  of  the  kernel  singularity,  is  outside  or  inside  the  boundary, 
respectively.  If  the  kernel  singularity  is  on  the  boundary  (Ct  €  dR),  a  is  equal  to  the  included 
angle  and  the  integral  is  treated  as  principal-valued.  For  most  of  our  computations,  we  take 
to  approach  the  boundary  from  outside  of  the  domain  so  that  a  is  zero.  The  algebraic  system 
is  formed  by  discretizing  the  integral  as  explained  in  section  3  and  letting  the  kernel  singularity 
approach  each  of  the  N  nodal  points,  C*  “*  A  special  limiting  process  is  then  needed  to 
evaluate  the  integration  near  Ct  • 

The  boundary  contour  dR  is  composed  of  dR4  and  dR^,  where  is  given  on  dR^  and 
is  given  on  dR^.  For  the  problem  to  be  well-posed,  either  ^  or  ^  must  be  specified  along 
the  entire  boundary.  This  is  related  to  the  Dirichlet  or  Neumann  problem,  respectively,  in 
the  standard  formulation.  Although  generalized  (Robin)  boundary  conditions  can  be  included 
in  these  computational  schemes,  we  do  not  study  them  here.  Hence,  the  boundary  can  be 
categorized  as  the  following  two  types: 

Re{0}  =  d  is  given  on  dR^  (3) 
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/m{/3}  =  w  is  given  on  dli^  . 


(4) 


When  the  boundary  condition  is  given  exclusively  as  either  0  oi  ii>  along  the  entire  boundary, 
the  solution  is  not  unique,  since  an  additional  real  or  imaginary  constant  can  be  added  to  the 
solution.  This  is  related  to  a  solvability  condition  for  the  Neumann  problem.  Therefore,  an 
appropriate  additional  condition  or  a  replacement  condition  (if  the  problem  is  not  to  be  overde¬ 
termined)  is  required.  We  have  found  that  replacing  the  boundary  condition  at  one  point  (i.e., 
assuming  <fi  is  known  at  every  point  except  one,  where  v>  is  known)  yields  much  less  accurate 
results.  Therefore,  we  use  an  additional  condition,  and  hence,  all  algebraic  systems  are  overde¬ 
termined,  although  sometimes  only  by  one  extra  equation. 


To  examine  the  corresponding  singularity  system  in  the  standard  formulation,  the  integrand 
in  (2)  can  be  rewritten  as 


ds  +  i 


v)e 


*-C» 


ds  , 


where  d  is  the  tangent  angle  of  the  boundary  with  the  positive  x  axis.  Hence,  the  contribution 
from  the  above  differential  element  to  the  complex  potential  inside  the  domain  is  given  by 

d>e<* 


•-Ck 


ds 


(5) 


This  relation  shows  that  the  integral  equation  (2)  is  equivalent  to  a  system  of  distributed 
normal  dipoles  of  strength  ^  and  tangential  dipoles  of  strength  on  the  boundtiry.  However, 
we  can  represent  the  potential  Reid  by  using  only  one  of  the  dipole  systems,  following  the  con¬ 
cepts  of  the  indirect  boundary  integral  method.  (In  the  indirect  method  [16],  the  computational 
unknowns  are  the  strengths  of  the  distributed  singularities  rather  than  the  original  solution,  as 
in  the  direct  method.)  Baker  et  al.  [8]  shows  that  dipole  distributions  are  more  stable  than 
source  or  vortex  distributions  for  nonlinear  breaking  wave  simulations.  We  could  formulate  a 
complex  integral  method  similar  to  a  combination  of  the  sources  and  vortices  by  introducing,  for 
example,  a  logarithmic  singularity  into  (2)  rather  than  the  simple  pole.  However,  our  preliminary 
numerical  studies  indicate  that  the  simple  pole  in  (2)  performs  better  than  other  singularities. 
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including  those  of  higher-order. 


3  Numerical  Model 

Following  the  procedure  of  Vinje  and  Brevig  [4],  both  the  boundary  contour  and  the  complex 
potential  are  represented  as  piecewise-linear  functions.  The  location  and  complex  potential  at 
the  j**  node  point  are  given  by  r,  and  0j,  respectively.  Then  0  on  panel  is  described  by 

0  —  T.  7*  “  ^)  +  -  *i)]  .  (6) 

‘i+i  “  V 

We  can  choose  C*  of  Eq.  (2)  to  be  anywhere,  although  since  we  treat  the  case  o  =  0.  it 
must  be  outside  the  solution  domain.  Typically,  it  is  chosen  to  approach  some  point  on  the  ik‘^ 
panel.  The  panel  midpoint  and  the  node  point  are  favored  equally  for  Ct  iii  previous  formulations. 
Our  studies  show  no  significant  differences  in  accuracy  and  convergence  rates  between  those  two 
cases.  Hence,  we  choose  the  endpoints,  which  are  easier  to  code.  We  also  have  studied  a  more 
overconstrained  problem  by  placing  control  points  at  the  node  points  and  the  midpoints  and  have 
found  no  significant  improvement.  In  section  9  we  examine  the  advantages  of  moving  sway 
from  the  contour. 


Letting  approach  each  node  from  the  outside  yields  the  algebraic  equivalent  of  (2)  approx¬ 
imated  as  [1] 


N 


=  0 


for  t  =  1, . . . ,  ^ 


(7) 


where 


filizi*  In  ( la  (  ^ 

*f+»  -  \  J  Zj-  *;-l  V  -  C*  / 


for  j  /  i. 


(8) 

(9) 


Eq.  (8)  is  evaluated  using  L’Hospital’s  rule  when  j  =  k  +  loik-l.  Moving  the  unknown 
boundary  conditions  to  the  right-hand  side  gives  the  following  complex  algebraic  system  for 
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unknown  6j  or 


53  +  i  53  =  -  53  13  ^jk4>j,  iot  k=l,...,N.  (10) 

iidRt  jiOR*  j€dR*  j€dR^ 

r^t  represents  the  influence  of  the  i:‘*  node  kernel  singularity  on  the  j"*  node  complex  potential 

3j.  In  general,  Ft-i,*,  Ftt,  and  Tk+i,*  have  the  dominant  influence  on  the  /b‘*  equation  except 

for  contours  that  are  not  convex  or  have  sharp  cusps.  The  imaginary  part  of  is  nearly  — x 

and  the  real  part  is  nearly  zero  for  uniformly-spaced  nodes  on  smooth  parts  of  the  contour,  and 

exactly  -r  and  zero,  respectively,  for  flat  contours.  On  the  other  hand,  the  imaginary  parts  of 

Fk_i,k  and  Ft+i,*  are  nearly  zero  while  the  real  parts  are  not  small. 

Choosing  the  imaginary  part  of  system  (8)  for  the  ifc‘*  equation  where  v»  is  given  and  the 
real  part  where  o  is  given,  results  in  a  system  matrix  with  a  strong  diagonal.  Taking  the  other 
systems  results  in  a  system  matrix  with  a  weak  diagonal.  We  now  refer  to  the  former  system 
as  the  strong  system  and  the  latter  as  the  weak  system.  The  strong  system  b  very  similar  to 
a  Fredholm  integral  equation  of  the  second  kind,  and  the  weak  system  is  similar  to  a  Fredholm 
integral  equation  of  the  first  kind.  As  indicated  in  the  introduction,  we  solve  the  weak  and  strong 
systems  together  as  an  overdetermined  system  in  a  least-squares  sense  and  compare  these  results 
to  those  obtained  from  the  individual  systems. 


4  Truncation  Error  in  Numerical  Integration 


The  truncation  error  in  the  numerical  integration  is  analyzed  for  the  case  of  uniformiy-sp2u:ed 
nodes.  In  general,  the  dominant  truncation  error  in  (2)  for  the  k'^  control  point  occurs  in  the 
integration  from  z  =  zk~i  to  z  =  Zk+i,  as  shown  in  Figure  1.  Rewriting  the  integral  (2)  as 


Jtu~l  *  ~  Ct 


(11) 


For  uniformly-spaced  nodes  with  slowly  varying  geometric  curvature,  we  can  examine  the  integral 
Ik  by  expanding  0{z)  at  z  =  (k  as, 


0(2)  =  0(Ck) + 0'«k)(z  -  ct) + -  a)' + o(iiz  -  uii^) .  (12) 


Using  (12),  we  can  write  Eq.  (11)  as 


Ik  =  ll^^  +  I^^^  +  li^^  +  0(\\z-Ckf)  , 


(13) 


where 


(14) 

/i‘^  =  ^'(C*)(^*+i-^t-i).  (15) 


and 

=  \0"«k)[(zk^i  -  n)*  -  (Zk-x  -  Zk)^  .  (16) 

Since  our  numerical  model  assumes  0  varies  linearly  between  nodes,  the  dominant  truncation 
error  in  the  calculation  of  /»  is 

When  0  is  linear,  the  integration  is  exact  and  the  integral  equation  algorithm  gives  exact  so¬ 
lutions.  In  Figure  2,  a  very  strong  correlation  between  /Ze{4’^}  and  the  solution  error  is  shown 
fw  a  problem  inside  a  circular  contour  with  uniformly-spaced  nodes  when  ^  is  given  as  /Ie{sin  z). 
The  correlations  are  nearly  1.0  and  -I.O  for  the  weak  and  strong  smems,  respectively.  Figure 
3a  shows  that  the  strong  and  weak  systems  give  approximately  the  same  errors  but  of  opposite 
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sign.  The  error  distribution  along  the  same  contour  for  another  lest  function,  in  which  3  —  e“”, 
is  shown  in  Figure  3b.  The  results  for  a  problem  when  the  solution  domain  is  outside  the  unit 
circle  are  given  in  Figure  3c.  This  example  represents  a  dipole  in  a  uniform  potential  gradient. 
The  imaginary  part  of  the  potential  should  be  zero  on  the  unit  circle.  To  solve  this  problem, 
the  uniform  part  must  be  subtracted  out  so  that  there  is  no  contribution  from  the  contour  at 
infinity.  A  branch  cut  is  chosen  to  connect  the  contour  at  infinity  to  the  inside  circular  contour, 
which  is  then  integrated  in  a  clockwise  sense.  All  examples  in  Figure  3  show  the  typical  error 
cancellation  characteristics  of  the  overdetermined  system.  The  overdetermined  system  averages 
the  errors  of  weak  and  strong  systems  so  that  the  error  distribution  curves  of  the  overdetermined 
system  nearly  collapse  to  the  abscissa.  Other  calculations  show  that  this  error  cancellation  for 
the  outer  problem  is  nearly  complete  when  the  dipole  is  moved  from  the  origin  but  still  well 
inside  the  contour. 

It  is  noteworthy  that  the  solutions  of  the  weak  system  sometimes  show  small  oscillations 
with  wavelengths  twice  the  nodal  spacing.  This  common  zig~zag  instability  is  due,  in  part,  to 
the  nearly  singular  matrix.  At  times,  when  this  numerical  instability  becomes  large,  the  weak 
system  matrix  becomes  so  numerically  singular  that  direct  inversion  is  no  longer  possible.  Then 
iterative  refinement,  as  done  here,  is  required. 

For  noncircular  contours,  the  errors  contributed  from  ofT-diagonal  elements  grow  as  the  local 
geometric  curvature  increases.  Figure  4  shows  the  error  distribution  curve  and  the  correlation 
between  error  and  /2e{/^^^}  when  9  is  given  for  an  dlipticai contour  with  an  eccentricity  of  0.91b5 
(the  major  axis  of  unit  radius  for  the  examples  shown  here  is  always  centered  on  the  x-axis).  A 
relatively  strong  correlation  still  exists,  but  the  error  of  the  weak  system  is  more  sensitive  near 
zero  In  this  case,  the  cancellation  characteristic  of  the  overdetermined  system  is  no  longer 
nearly  perfect,  as  shown  in  Figure  3.  Table  I  shows  the  correlation  numhers  between  the  solution 
error  and  for  various  ellipses  with  eccentricities  of  0.0  to  0.9798  (minor  axes  in  the  y 


direction  of  1  to  0.1).  The  weak  system  loses  its  strong  correlation  at  an  eccentricity  of  0.9798, 
where  large,  spurious  oscillatory  behavior  develops  in  the  solution. 

We  will  show  that  the  error  cancellation  of  the  overdetermined  system  seems  more  dependent 
on  the  variation  of  geometric  curvature  and  on  the  nodal  spacing,  rather  than  on  the  solution 
characteristics. 
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5  Numerical  Investigation 


Comparisons  of  accuracy,  convergence  rate,  and  computing  time  are  made  for  weak,  strong, 
and  overdetermined  sets  of  system  ( 10).  The  effects  of  nonuniform  nodal  spacing,  types  of  bound¬ 
ary  conditions,  nearly  singular  solution  characteristics,  and  contour  shape  on  the  solution  errors 
are  also  investigated. 

The  overdetermined  matrix  problem  can  be  solved  two  ways.  One  method  uses  a  routine 
based  on  the  Householders  transformation  [9].  This  routine  takes  an  order  of  V  operations, 
where  U  is  the  number  of  unknowns  and  V  is  the  number  of  equations.  The  other  uses  a  conju¬ 
gate  gradient  iterative  technique  [10],  which  solves  the  problem  in  an  order  of  UV'L  operations, 
where  L  is  the  number  of  iterations.  The  computational  savings  of  the  iterative  technique  would 
be  important  if  time  marching  were  desired,  especially  since  a  good  initial  guess  is  available  from 
the  previous  time  step.  However,  even  with  a  homogeneous  initial  guess,  a  typical  solution  for  a 
nonsingular  matrix  requires  less  than  10  iterations. 

The  iterative  solver  often  takes  more  computing  time  than  the  Householders  transformation 
method  for  the  weak  system,  but  iterative  methods  are  more  reliable  when  the  system  is  nearly 
singular,  as  is  true  in  the  weak  system.  Here,  we  always  use  the  conjugate  gradient  iterative 
method  with  a  zero  initial  iterate  to  compare  the  efficiency  of  the  three  systems. 

5.1  Rate  of  Convergence  and  Computing  Time 

The  root  mean  square  errors  (£])  and  corresponding  computing  times  are  examined  as  a 
function  of  the  number  of  nodes  to  determine  convergence  rates  and  computational  efficiency. 
The  maximum  errors  (£00)  follow  a  very  similar  pattern  and  hence  are  not  shown.  Figure  5 
shows  £2  for  the  three  systems  in  log-log  scale  with  uniformly-spaced  nodes  on  a  circular  contour 
when  ^  is  given  as  /fe{8inz}.  Single-precision  (six  digits)  results  are  limited  to  approximately 
five  decimal  point  accuracy  as  round-off  error  dominates  truncation  error  beyond  N=32.  The 
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overdetermined  system  is  more  accurate  by  one  digit  at  ^  =  16,  and  the  convergence  rate  of  the 
overdetermined  system  is  twice  that  of  the  weak  or  strong  systems. 

In  Figure  6,  computing  times  of  the  single^precision  program  are  shown.  For  the  same  N,  the 
computing  time  of  the  overdetermined  system  is  the  largest  and  the  computing  time  to  solve  the 
weak  system  increases  rapidly  when  >  32  as  the  system  matrix  becomes  more  singular  and 
the  number  of  iterations  increases.  Figure  6  clearly  shows  that  the  overdetermined  solution  is 
the  most  efficient  of  the  three  systems  for  this  example. 

The  effect  of  the  contour  shape  on  the  efficiency  of  inner  problems  is  examined  for  elliptical 
contours  with  various  eccentricities  (e  =  0.0  to  0.995)  for  the  same  0.  The  diameter  of  the  major 
axis  is  fixed  to  be  2.0  and  centered  on  the  x  axis.  .\s  the  eccentricity  increases,  the  Ej  of  the 
overdetermined  system  increase,  but  those  of  the  strong  systems  decrease,  as  shown  in  Figure 
7.  The  £]  of  the  weak  system  for  «  >  0.9165  increase  due  to  the  more  oscillatory  behavior  of 
the  solution.  In  Table  II,  the  £2  are  given  as  a  function  of  the  number  of  nodes  for  various  el¬ 
lipses.  In  most  cases,  the  overdetermined  system  gives  the  best  results,  except  for  highly  eccentric 
elliptical  contours  where  the  weak  system’s  oscillations  affect  the  overdetermined  solution  as  well. 

The  £2  convergence  rates  determined  at  the  largest  computed  N  are  compared  in  Table  II 
for  inner  problems  for  ellipses  with  various  eccentricity.  Astonishingly,  the  convergence  rate  of 
the  overdetermined  system  for  a  circular  contour  (e  =  0)  approaches  4.0.  The  convergence  rates 
of  all  the  other  examples  are  the  predicted  2.0. 

£2  for  outer  problems,  when  ^  s  1/z,  are  also  calculated  and  given  in  Table  III  for  various 
ellipses  (0  <  e  <  0.995)  with  a  double-precision  program.  The  results  of  the  overdetermined  sys¬ 
tem  are  the  best  for  all  cases  tested.  The  convergence  rates  approach  2.0  for  the  weak  and  strong 
systems.  The  convergence  rates  of  the  overdetermined  system  also  approach  2.0  except  at  e  =  0, 
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which  again  exhibits  fourth-order  convergence.  We  also  have  obtained  fourth-order  accuracy  for 
exterior  problems  when  the  simple  poles  are  placed  in  the  interior,  not  close  to  the  boundary. 
As  the  eccentricity  increases,  Ej  increases  for  all  systems  due  to  the  more  rapid  variation  of  the 
function  and  the  kernel. 

Since  our  numerical  model  uses  a  piecewise-linear  representation  for  the  potential  function,  we 
expect  second-order  convergence.  However,  for  some  cases,  the  convergence  rates  of  the  overde¬ 
termined  system  approach  fourth-order  due  to  the  cancellation  of  errors  between  the  weak  and 
the  strong  systems.  Our  examples  show  that  these  cancellations  occur  when  the  contour  is  circu¬ 
lar  (or  nearly  so)  with  uniformly-spaced  nodes.  Some  calculations  not  described  here  have  shown 
third-  and  sixth-order  convergence  for  the  strong  and  overdetermined  systems,  respectively,  when 
a  cubic-spline  representation  is  used  for  i3. 


5.2  Nonuniform  Nodal  Spacing 

Sometimes  nonuniform  nodal  spacing  is  desirable  when  solution  refinement  is  needed  due  to 
large  solution  gradients.  Also,  the  nodal  spacing  in  problems  with  convective  nodes  (such  as  free 
surface  flow  problems)  becomes  nonuniform.  This  can  cause  numerical  instabilities,  so  that  a 
Altering  scheme  or  a  regridding  scheme  is  needed  to  maintain  uniform  spacing  [11]. 

In  Figure  8,  a  schematic  of  nonuniformly-spaced  nodes  is  shown.  The  nonuniformity  at  node 
k  (in  an  otherwise  uniformly-spaced  grid)  is  defined  as  the  ratio  between  L,  and  Lq.  In  this  case, 
the  real  part  of  Ftt  is  no  longer  small.  Moreover,  the  dominant  off-diagonal  terms  and 

rt4.i,i  also  change  significantly  while  the  imaginary  part  of  Tkk  remains  around  —v.  Hence,  the 
local  nonuniformity  greatly  affects  the  global  solution  characteristics. 

Figure  9  shows  the  same  example  as  in  Figure  3a.  except  that  node  1  is  moved  to  form  various 
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nonuniformities  as  defined  by  Figure  8.  The  errors  of  the  strong  system  are  affected  only  near 
that  nonuniformly-spaced  node.  As  the  nonuniformity  increases,  the  solution  of  the  weak  system 
deteriorates  globally  due  to  the  zig-zag  instability.  However,  if  a  numerical  filter  were  used  to 
eliminate  oscillations  with  wavelengths  equal  to  the  node  spacing,  the  weak  system  error  is  then 
comparable  with  the  strong  system  and  again  of  opposite  sign.  This  is  why  the  solution  for 
the  overdetermined  system  remains  very  good.  Note  that  the  error  curve  of  the  overdetermined 
system  cannot  be  distinguished  from  the  zero  axis  except  near  the  moved  node. 

The  correlation  between  and  the  solution  error  is  still  strong  for  strong  systems,  except 
around  the  node,  but  it  is  very  poor  for  the  weak  system.  The  computing  time  for  the  weak 
system  increases  rapidly  as  the  nonuniformity  increases,  since  the  number  of  iterations  grows 
rapidly  as  the  matrix  becomes  more  numerically  singular. 

5.3  Types  of  Boundary  Condition 

When  the  boundary  contour  b  composed  of  both  dRf  and  diZ^,  the  intersecting  node  can 
present  a  difficulty  even  without  any  geometric  discontinuity.  We  could  expect  the  solution  to 
be  singular  at  the  intersecting  node.  We  also  would  usually  expect  that  the  change  in  boundary 
conditions  to  occur  at  a  comer  in  the  boundary  contour  (which  brings  other  difficulties  that  will 
be  discussed  in  section  5.5).  We  do  not  consider  either  case  here. 

We  examine  a  circular  contour  with  ^  given  on  the  upper  half  (y  >  0)  and  0  on  the  lower 
half.  Figure  10  shows  the  error  distribution  curves  for  two  cases  when  0  =  e~** .  Either  b  or  tf*, 
or  both  can  be  treated  as  known  at  the  intersecting  nodes  (y  =  0).  Only  one  boundary  condition 
is  prescribed  at  each  node  in  (A),  while  both  i  and  d*  are  given  at  the  intersecting  nodes  in  (B). 
In  both  cases,  the  error  of  the  overdetermined  system  is  much  smaller  while  the  weak  system 
solution  again  shows  oscillatory  behavior. 


£2  for  these  two  cases  are  compared  in  Figure  11.  The  solutions  when  the  additional  known 
boundary  constraints  are  used  are  better  for  ail  three  systems.  Convergence  rates  remain  the 
same  and  the  improvement  in  accuracy  is  small.  There  is  also  no  large  local  error  at  the  inter¬ 
section  point.  We  conclude  that  the  change  in  boundary  condition  type  by  itself  is  not  difficult 
to  handle  computationally. 

5.4  Singularity  Near  Contour 

When  the  boundary  condition  values  change  abruptly  due  to  nearby  singularities,  the  solution 
error  at  nodes  near  those  singularities  is  expected  to  increase.  As  an  example,  the  error  for  a 
unit  radius  circular  contour  where  the  solution  has  a  simple  pole  near  the  contour  is  examined, 
i.e,  o  is  given  as  the  real  part  of 

0= - - - 

z-(l  +  0’ 

where  e  is  a  small  constant.  The  solution  errors  for  rjf  of  the  strong  system  increase  mainly  at 
those  nearby  nodes,  whUe  those  of  the  weak  system  exhibit  an  oscillatory  behavior  with  peaks  at 
nearby  nodes,  as  shown  in  Figure  12.  As  e  increases,  the  oscillation  aunplitude  diminishes  and  the 
weak  system  finally  recovers  the  same  error  pattern  as  the  strong  system,  but  again  of  opposite 
sign. 

E}  shown  in  Figure  13  show  the  advantage  of  the  overdetermined  system  in  accuracy  and 
convergence  rate.  The  convergence  rate  for  the  overdetermined  system  is  again  4.0  and  those 
of  the  other  systems  are  around  2.0.  For  this  function,  the  nodal  spacing  must  be  less  than  e 
before  the  high  convergence  rate  is  achieved.  We  note  that  the  weak  system  is  unstable  due  to  an 
abrupt  change  in  which  again  implies  that  the  system  matrix  of  the  weak  system  is  nearly 
singular. 
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5.5  Contours  with  Corners 

The  boundary  integral  accuracy  is  known  to  deteriorate  in  the  vicinity  of  contour  corners 
due  to  high  geometric  curvature.  Moreover,  the  solution  is  usually  singular  at  corners.  Grisvard 
[12]  has  the  most  mathematically  complete  discussion  of  numerical  singularities  arising  from  the 
nonsmooth  domain  for  elliptic  boundary  value  problems.  Pina,  Fernandes,  and  Brebbia  [13] 
suggest  using  mesh  refinements  to  overcome  the  difficulties  of  nonoptimal  convergence  arising 
from  an  irregular  boundary  contour.  These  refinements  become  more  important  if  the  solution  is 
known  to  be  singular  or  to  vary  rapidly  at  a  corner  (Lin,  Newman  and  Yue  [14]).  The  numerical 
model  can  be  truly  singular  if  a  modified  variable  is  formed  by  subtracting  out  the  local  singular 
solutions  from  the  original  variable  (Kelmanson  [7]).  Mesh  refinements  near  corner  can  create 
difficulties  with  numerical  instability  such  as  at  the  free  surface  of  the  wave-maker  problem  [11]. 

Even  if  the  solution  is  not  singular  at  a  corner,  error  is  induced  in  large  part  by  the  contri¬ 
bution  from  the  abrupt  change  in  the  geometric  curvature.  This  can  be  seen  in  Eq.  (16)  for 
At  a  corner,  the  term  (zt+x  “  C*)*  —  (^s-i  —  C*)*  >*  not  small.  If  0"{^k)  is  nonzero  at  the 
comer,  the  error  in  numerical  integration  will  be  much  larger  there.  In  this  case,  the  error  has 
its  maximum  near  the  corner  and  can  spread  to  the  entire  contour. 

A  modified  circular  contour  with  a  corner  (Figure  14)  is  used  to  test  the  effects  of  the  corner 
on  the  solution  accuracy.  In  this  case,  the  imaginary  part  of  Tkk  at  the  cmner  stays  at  — t/2 
as  the  number  of  nodes  increases.  However,  the  off-diagonal  influence  coefficients  are  strongly 
affected  by  the  comer  and  any  refinement  there. 

As  Figure  15a  shows,  the  error  of  the  strong  system  is  concentrated  at  the  corner  while  it 
remains  at  the  previoxis  error  leveb  on  the  smooth  part  of  the  contour,  as  shown  in  Figure  3a. 
The  error  of  the  weak  system  solution  oscillates  over  the  entire  contour,  and  the  computing  time 
(the  nundxer  of  iterations)  increases  rapidly  due  to  the  corner.  However,  if  a  numerical  filter 
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were  used  to  eliminate  oscillations  with  wavelengths  equal  to  the  node  spacing,  the  weak  system 
would  actually  have  smaller  error  than  the  strong  system.  The  solution  of  the  overdetermined 
system  is  still  significantly  better,  but  the  maximum  error  is  one  thousand  times  greater  than 
the  solution  of  the  overdetermined  system  for  the  smooth  contour,  shown  in  Figure  3a.  All  three 
approaches  show  second-order  convergence  for  contours  with  corners. 

The  given  solution  for  this  example  is  analytic  at  the  corner,  but  we  still  have  very  grave 
computational  difficulties  there.  To  make  smaller,  we  develop  a  modified  complex  potential 
defined  by 

(17) 

where  is  the  coordinate  of  the  corner.  Thus,  the  second  derivative  of  the  modified  complex 
potential  dc  is  forced  to  be  zero  at  the  corner.  This  eliminates  the  highest-order  error  term  in 
our  piecewise-linear  integration  routine  (Eq.  16)  at  the  corner. 

Figure  15b  shows  improved  results  using  this  nx>dification,  especially  at  the  corner.  All  six 
solutions  (we^,  strong,  and  overdetermined  systems  for  both  the  modified  and  unmodified  for¬ 
mulations)  still  appear  to  be  second-order  convergent  for  contours  with  corners.  However,  the 
£3  or  £00  for  fixed  N  is  significantly  improved  by  the  modified  method.  For  this  example,  with 
iV  =  64,  these  errors  are  reduced  by  a  factor  of  50. 

These  results  suggest  using  an  iterative  method  for  more  accurate  solutions  to  nontesi  prob¬ 
lems,  when  the  unknown  value  of  0"(,ie)  i*  updated  to  a  more  correct  value.  Further  developments 
of  this  method  are  required  for  problems  having  more  than  one  corner  and  for  the  outer  problems. 
It  is  also  expected  that  even  more  accurate  solutions  can  be  obtained  by  requiring  the  higher 
derivatives  of  the  modified  complex  potential  to  be  zero  at  corners. 
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6  Desingularization 

Patterson  and  Sheikh  [15]  introduce  a  regular  boundary  integral  method  by  moving  the  sin¬ 
gularity  of  the  fundamental  solution  away  from  the  boundary  contour.  Their  rationale  for  doing 
this  appears  to  be  chiefly  for  desingularizing  the  kernel.  Hence,  all  integrals  can  be  computed 
numerically.  They  also  claim  the  method  “tolerates  higher-order  singularities”  and  base  these 
conclusions  on  comparisons  with  finite  element  methods  for  elastostatic  problems. 

We  show  here  that  moving  Xi  in  Eq.  (2)  well  outside  the  boundary  has  several  advantages 
over  the  standard  approach.  While  we  do  not  perform  the  integrals  by  simple  quadrature  (we  still 
use  Eq.  (fl)),  it  is  clear  that  numerical  integration  would  no  longer  seriously  affect  the  results. 
The  determination  of  the  influence  coefficients  is  still  simpler,  since  special  cases  are  not  required 
(as  when  integration  passes  through  a  singular  point  of  the  kernel).  Now.  however,  we  must 
decide  where  to  locate  the  kernel  singularities. 

We  use  two  strategies  for  locating  the  kernel  singularities  outside  the  contour  shown  in  Figure 
16.  The  first  strategy  uses  the  perpendicular  bisector  of  the  straight  contour  between  adjacent 
nodes.  This  strategy  has  two  potential  difficulties.  If  f  is  too  small,  the  singular  point  may  lie 
inside  an  actual  convex-curved  contour.  If  /  is  too  large,  the  singular  point  may  lie  inside  the 
domain  of  a  highly  contorted  contour.  The  second  strategy  places  all  singular  points  on  a  circle 
whose  center  is  roughly  at  the  centroid  of  the  contour.  This  strategy  avoids  the  difficulties  of 
the  first,  but  it  may  not  place  the  kernel  singularities  sufficiently  close  to  the  contour  to  capture 
nearly  singular  behavior  of  the  solution.  The  second  strategy  also  has  the  freedom  of  placing  any 
number  of  kernel  singularities  on  the  circle,  and  hence,  the  algebraic  system  can  become  more  or 
less  overdetermined. 

Figure  17  shows  the  error  computed  using  the  first  strategy  and  two  values  of  /;  in  this 
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case  the  nodes  are  placed  evenly  on  a  circular  contour  with  0  =  sin  z.  Here,  we  see  that  the 
convergence  is  exponential  and  that  the  overdetermined  system  is  again  the  best.  Figure  17 
indicates  that  placing  the  singularities  further  from  the  contour  is  preferable.  Figure  18  shows 
this  to  be  the  case  for  the  same  example  with  fixed  N  =  16.  However,  a  price  is  ptud  in  that 
the  algebraic  system  becomes  less  diagonally  dominant  as  /  increases,  requiring  more  iterations 
and  a  tighter  tolerance  on  the  convergence  parameter  to  achieve  the  acquired  accuracy.  For  this 
contour  we  note  that  when  /  <  (1  —  cos  ^)/ sin  y  as  0.05  for  N  =  16,  the  singular  point  ac¬ 
tually  moves  inside  the  contour.  Figure  18  shows  this  happens  without  greatly  changing  the  error. 

Figure  19  also  shows  exponential  convergence  for  the  circular  contour  when  there  is  a  solution 
singularity  outside  the  domain.  In  this  case,  a  semi-log  plot  indicates  the  true  exponential  char¬ 
acter  of  the  convergence.  (Only  the  overdetermined  results  are  shown.)  The  nodal  spacing  must 
be  small  compared  to  c  for  exponential  convergence  to  be  achieved.  As  a  result,  the  convergence 
when  e  =  0.1  (not  shown)  is  very  poor  ail  the  way  to  N  =  64.  As  might  be  expected,  the  most 
accurate  solutions  are  no  longer  obtained  when  /  is  large,  since  placing  the  kernel  singularities 
further  away  cannot  capture  the  effect  of  solution  singularities  close  to  the  contour. 

Although  the  overdetermined  system  still  has  smaller  errors,  it  is  not  significantly  smaller 
than  the  errors  for  the  desingularized  version  of  the  strong  or  weak  system  (Figure  13).  Simi¬ 
larly,  Figure  20  shows  that  the  same  is  true  for  elliptical  contours.  As  before,  all  systems  have 
second-order  convergence,  but  the  overdeternnined  system  is  slightly  more  accurate  than  the 
strong  or  weak  system.  (The  weak  systera  Is-not  shosro-since  its  results-are-very  similar  to  the  - 
strong  system.)  While  desingularization  greatly  improves  the  accuracy  of  the  weak  and  strong 
system,  it  only  slightly  improves  the  accuracy  of  the  overdetermined  system. 

For  the  examples  in  Figures  17-20,  the  desingularized  method  is  better  than  the  method  of 
Vinje  and  Brevig.  However,  both  strategies  for  placing  the  kernel  singularities  off  the  contour  do 
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not  yield  as  accurate  results  when  the  contour  has  a  corner.  The  convergence  rate  is  still  the  same 
and  the  absolute  errors  increase  slightly  (not  shown,  but  typically  25  percent).  Further  work  is 
required  to  determine  how  to  place  the  singular  points  next  to  the  contour  corners;  however, 
our  computations  show  that  decreasing  /  or  R,  near  the  corners  is  not  beneficial,  even  when  no 
solution  singularity  exists  there. 

The  second  strategy  shown  in  Figure  16  allows  an  arbitrary  number  of  kernel  singularities  to 
be  used  in  formulating  the  algebraic  system.  When  the  real  and  imaginary  systems  are  used,  the 
number  of  singular  points,  N,,  can  be  less  than  the  number  of  nodes.  JV,  actually  N,  >  iV/2. 
Table  IV  shows  typical  error  results  using  variable  N,.  All  results  are  for  a  unit  circle  contour. 
.V  =  16,  and  0  =  sinz,  with  N,  singular  points  equally  spaced  on  R,  =  1.5.  Table  IV  shows 
that  S,  can  be  reduced  below  V;  increasing  S,  further  increases  the  accuracy  insignificantiy. 
However,  increasing  N,  for  the  strong  equations  to  overconstrain  the  system  can  be  beneficial  if 
Nt  is  twice  N.  Table  IV  does  not  show  that  overdetermining  the  algebraic  system  in  this  manner 
is  more  expensive;  more  logarithms  must  be  computed  to  form  the  influence  the  matrix  and  the 
resulting  system  is  not  as  diagonally  dominant  (i.e.,  more  iterations  are  required  to  solve  the 
algebraic  system). 


7  Spectral  methods 


Following  the  method  of  Baker  (17],  exponential  convergence  can  also  be  achieved  by  inte¬ 
grating  the  principal-valued  form  of  (2)  using  the  trapezoidal  rule.  To  obtain  this  convergence, 
the  integral  must  be  periodic  and  integrated  with  respect  to  a  real  variable,  in  this  case,  an 
arciength  parameter.  The  singular  part  of  the  integral  is  avoided  by  evaluating  the  integrand  at 
every  other  point.  The  integral  then  becomes 


/ 


d{$)  dx 
r(s)  -  <;*  d$ 


d$  =  iaffiCk)  . 


(18) 
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Baker  indicates  this  expression  will  give  exponentially  convergent  results  if  the  shape  (and 
hence  dz/ds)  is  known  exactly.  However,  in  contour  dynamics  problems,  for  example,  the  shape 
is  not  known  precisely.  This  problem  can  be  overcome  if  the  contour  is  known  to  be  smooth,  since 
c(s)  and  hence  dzfds  can  then  be  determined  with  exponential  accuracy  using  spectral  methods. 
We  have  performed  these  calculations  using  a  standard  fast  Fourier  transform  routine  (FFT)  for 
both  z(s)  and  j/(s)  to  obtain  exponential  accuracy  for  fully  nonlinear  gravity  waves.  These  results 
will  be  reported  elsewhere.  Here,  we  present  results  for  circles,  ellipses,  and  squares,  where  the 
shape  is  represented  by  the  same  Fourier  representation.  For  the  elliptical  contours,  the  nodes 
are  placed  at  constant  polar  angle  intervals  as  measured  from  the  center  of  the  ellipse.  Then  the 
FFT  is  exact  and  the  convergence  is  exponential,  with  the  highest  convergence  for  nearly  circular 
contours,  as  shown  in  Figure  21  for  i3  =  sinz.  The  overdetermined  system  results  are  nearly  the 
same  as  the  other  systems  and  hence  only  one  is  shown. 

A  simple  problem  where  the  FFT  is  not  exact  is  given  using  the  same  example  of  the  ellipse 
(e  =  0.866)  with  nonuniform  nodal  spacing.  Figure  22  shows  the  convergence  for  the  case  when 
the  node  is  moved  by  0  =  .002T/Ar.  Even  for  this  nearly  imperceptible  movement,  the 
asymptotic  convergence  for  all  three  systems  is  first  order.  The  overdetermined  system  yields 
slightly  better  Ei-  The  error  distribution  curves  for  the  three  systems  are  shown  in  Figure  23 
for  iV  =  32.  While  £}  is  nearly  the  same  for  all  systems,  the  maximum  error  (£«)  for  the 
overdetermined  system  is  significantly  lower.  All  three  systems  show  an  oscillatory  error  with 
the  wavelength  on  the  order  of  the  nodal  spacing.  This  could  give  rise  to  the  zig-zag  instability 
found  in  many  contour  dynamics  computations  [8].  The  oscillatory  error  of  the  overdetermined 
solution  looks  sufficiently  regular,  such  that  a  simple  filter  may  recover  exponential  convergence. 

Finally,  we  show  the  convergence  of  the  trapezoidal  method  in  Figure  24  for  a  square  contour 
centered  about  the  origin  with  each  side  of  length  two,  for  Q  =  sin  z.  We  see  in  this  case  that 
the  convergence  rate  is  much  less  than  first  order,  with  the  o\'erdetermined  system  giving  the 
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best  results.  Unlike  the  previous  example,  no  simple  filter  would  improve  the  results.  Part  of  the 
reason  for  the  poor  convergence  is  the  poor  convergence  of  d<^/ds  by  the  FFT  due  to  the  squares 
nonsmooth  contour.  Even  when  an  exact  d(/ds  is  used  (by  central  difference),  the  convergence 
(although  better)  is  still  less  than  first  order  due  to  the  discontinuities  in  d^/ds  at  the  corner 
nodes.  Any  hope  of  recovering  exponential  convergence  lies  in  moving  all  nodes  by  one-half  the 
nodal  spacing  to  avoid  placing  nodes  at  the  corners. 

It  should  be  noted  that  the  condition  of  the  matrix  obtained  in  this  manner  is  excellent,  and 
half  of  the  coefficients  are  zero.  This  means  that  the  algebratic  solver  saves  computational  effort, 
in  addition  to  the  time  saved  by  not  requiring  the  evaluation  of  logarithms.  The  additional  time 
to  compute  the  FFT's  is  not  significant. 

Obviously,  more  algorithm  development  is  needed  for  contours  with  corners.  In  addition,  it 
is  not  entirely  clear  how  to  extend  this  approach  to  domains  with  contours  that  are  not  simply 
connected. 

8  Conclusions 

Error  analysis  of  complex  variable  boundary  integral  methods  is  made  by  truncation  error 
analysis  and  numerical  investigation.  The  overdetermined  system  is  the  most  efficient  and  accu¬ 
rate  in  almost  all  cases;  this  sytem  averages  errors  of  the  opposite  sign  for  the  solutions  of  the 
strong  and  weak  systems.  The  overdetermined  system  has  a  fourth-order  convergence  rate  for 
circular  contours,  while  the  other  systems  exhibit  second-order  convergence.  When  the  contour  is 
not  circular,  the  convergence  for  all  three  systems  is  second  order,  but  the  overdetermined  system 
is  significantly  more  accurate  for  moderate  N.  All  systems  have  difficulties  at  contour  corners, 
which  can  be  alleviated  by  introducing  a  modified  complex  potential  such  that  the  derivatives 
are  zero  at  the  corner.  The  overdetermined  system  takes  the  fewest  number  of  iterations  in  our 
coqjugate  gradient  iteration  procedure.  As  a  result,  the  computational  time  is  less  than  twice 
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that  of  the  weak  or  strong  system.  The  increased  accuracy  of  this  method  more  than  repays  this 
additional  cost. 

For  smooth  contours,  exponential  convergence  can  be  achieved  by  desingularization  of  the 
kernel  or  by  using  a  trapezoidal  integration  rule.  The  exponential  convergence  is  destroyed  for 
noncircular  contours  when  desingularization  is  used  or  when  nodes  are  unevenly  placed  for  ei¬ 
ther  method.  The  spectral,  trapezoidal  method  is  much  more  accurate  if  the  nodes  are  precisely 
placed.  However,  small  errors  can  render  the  convergence  for  the  trapezoidal  method  to  be  first- 
order.  while  the  desingularization  method  is  more  robust.  Both  of  these  modifications  perform 
poorly  when  a  contour  corner  is  present.  The  overdetermined  system  is.  for  the  most  part,  the 
most  efficient  for  these  modifications  as  well. 

The  same  approach  may  be  extended  to  three-dimensional  potential  problems  by  using  two 
sets  of  integral  equations,  one  of  which  is  derived  from  Green's  theorem  and  the  other  from  the 
normal  derivative  of  Green’s  theorem.  Hence,  the  advantages  of  solving  an  overdetermined  sys¬ 
tem  may  yield  greater  efficiency  and  accuracy  for  threenlimensional  problems  as  well. 
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Figure  4:  Error  Distribution  and  Correlations  between  Error  and 
(ellipse.  e=0.9165,  inner  problem.  ^  =  sinr,  o  is  given 


Figure  5;  Ej  Errors  for  Various  Number  of  Elements 
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Figxire  10:  Error  Distribution  for  a  Mixed  type  Problem 
(circle, inner  problem,  /?  =  c~'*) 
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Figure  11:  £j  Errors  of  Mixed  Tj-pe  Problems 
(circle,  rs=l,  inner  problem,  0  =  c~'* 
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Figure  12:  Error  Distributions  for  Singular  Function  Case 
(circle,  inner  problem,  (p  is  given) 
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Figure  13:  Errors  for  Singular  Function  Case 
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Figure  14:  A  Modified  Circular  Contour  with  a  Corner 
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Figure  15;  Error  Distributions  for  a  Corner  Problem 
(inner  problem.|5  =  sin{z),  p  is  given) 
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Figure  18:  Effect  of  Desingiarization  Parameter  on  Error 


Figure  21:  Convergence  of  Trapezoidal  Method  for  nearly  Circular  Contours 
(ellipses,  inner  problem,  0  —  sin(z),  ^  is  given) 
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Figure  22:  Convergence  of  Trapezoidal  Method  for  Unevenly  Spaced  Node  Case 
(ellipse,  e  =  0.866,  0  =  sin(z)) 
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Figure  23:  Error  Distribution  of  Trapezoidal  Method  for  Unevenly 
Spaced  Node  Case 
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